# a function that runs RAS balancing procedure on a matrix A with vectors of constraints col and row
# NOTE: this version has pre-defined tolerance 0.0001 (to avoid manually setting tolerance in loops)
# ENSURE that Octave can find the function, use the "addpath" command: addpath ("/[insert the path to your Octave folder]/Octave")

# function arguments:
# name - initial estimate (matrix)
# col - column vector of contraints
# row - row vector of constrains

function A = ras (name, col, row)
# create a copy of the initial estimate that will be modifiied throughout the procedure
A0 = name;
do
# calculate a column vector of row multipliers
	RAS_r = col./sum(A0, 2);
# some elements in RAS_r may be NaN because of division of zero by zero, so replace NaN with 0
	RAS_r(isnan(RAS_r)==1) = 0;
# scale A0 row-wise
	A0 = diag(RAS_r)*A0;
# check the convergence column-wise
	res_s = row - sum(A0, 1);
# break the procedure if it converges
	if (max(abs(res_s)) <= 0.0001)
		A = A0;
		break;
	endif
# calculate a row vector of column multipliers
	RAS_s = row./sum(A0, 1);
# some elements in RAS_s may be NaN because of division of zero by zero, so replace NaN with 0
	RAS_s(isnan(RAS_s)==1) = 0;
# scale A0 column-wise
	A0 = A0*diag(RAS_s);
# check the convergence row-wise
	res_r = col - sum(A0, 2);
# condition
until (max(abs(res_r)) <= 0.0001)
A = A0;
endfunction
